HJB $\log K, R, Y$¶

\begin{align*} d (K_{d} + K_g) =& \left([\alpha_{d}+i_{d}-\frac{\phi_{d}}{2}i_{d}^{2}] K_d + [\alpha_{g}+i_{g}-\frac{\phi_{g}}{2}i_{g}^{2}] K_g \right) d t+\left(\sigma_{d}K_d + \sigma_g K_g \right) d W\\ d K / K =& \left([\alpha_{d}+i_{d}-\frac{\phi_{d}}{2}i_{d}^{2}] (1 - R) + [\alpha_{g}+i_{g}-\frac{\phi_{g}}{2}i_{g}^{2}] R \right) d t+\left(\sigma_{d}(1 - R)+ \sigma_g R \right) d W \\ d \log K =& \left([\alpha_{d}+i_{d}-\frac{\phi_{d}}{2}i_{d}^{2}] (1 - R) + [\alpha_{g}+i_{g}-\frac{\phi_{g}}{2}i_{g}^{2}] R - \frac{\mid\sigma_{d}(1 - R)+ \sigma_g R \mid ^2}{2} \right) d t\\ &+\left(\sigma_{d}(1 - R)+ \sigma_g R \right)d W \\ d R =& \left( (\alpha_g + i_g - \frac{\phi_{g}}{2}i_g^2) - (\alpha + i_d - \frac{\phi_d}{2} i_d^2)\right) \left[ R(1 - R)\right] dt + R^{2}(1-R)^{2}\left[\sigma_{d}^{2}+\sigma_{g}^{2}\right] d W \end{align*}\begin{align*} \delta v(\log K, R, Y; A_g') & =\max_{i_{g},i_{d} }\delta \log((A_d - i_d) (1 - R) + (A_g' - i_g)R ) + \delta \log K \\ & +\left(\{\alpha_{d}+i_{d}-\frac{\phi_{d}}{2}i_{d}^{2}\} (1 - R)+\{\alpha_{g}+i_{g}-\frac{\phi_{g}}{2}i_{g}^{2}\} R - \frac{\mid\sigma_{d}(1 - R)+ \sigma_g R \mid ^2}{2}\right) v_k \\ &+\frac{\mid \sigma_{d} (1 - R) + \sigma_g R\mid^2}{2}v_{kk}\\ &+ \left( (\alpha_g + i_g - \frac{\phi_{g}}{2}i_g^2) - (\alpha + i_d - \frac{\phi_d}{2} i_d^2)\right) \left[ R(1 - R)\right] v_R \\ &+\frac{1}{2}R^{2}(1-R)^{2}v_{RR}\left[\sigma_{d}^{2}+\sigma_{g}^{2}\right]\\ & +\beta_{f} (\eta A_d (1 - R) K )v_{Y}+\frac{1}{2}\varsigma^{2}(\eta A_d (1 - R) K)^{2}v{}_{YY} \\ &- [\{\gamma_{1}+\gamma_{2}Y_{t}\}\beta_{f}(\eta A_d(1 - R)K)+\frac{\gamma_2}{2}\varsigma^{2}(\eta A_d (1 - R) K)^{2}] \end{align*}

$\eta = 0.008$

In [1]:
# for post jump HJB
import os
import numpy as np
import pandas as pd
import pickle
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d
import plotly.graph_objects as go

mpl.rcParams["savefig.bbox"] = "tight"


with open("data/res-28-14-24", "rb") as file:
      data = pickle.load(file)

Kd = data["Kd"]
Kd_max = data["Kd_max"]
Kg = data["Kg"]
Kg_max = data["Kg_max"]
Kd_mat = data["Kd_mat"]
Kd_max = data["Kd_max"]
Kg_mat = data["Kg_mat"]
Kg_max = data["Kg_max"]
Y = data["Y"]
Y_mat = data["Y_mat"]
i_d = data["i_d"]
i_g = data["i_g"]
v = data["v0"]

print(v.shape)

(36, 99, 41)
In [2]:
fig = go.Figure(data=[go.Surface(z=v[:,:,0], x=np.exp(Kd_mat[:,:,0]), y=Kg_mat[:,:,0] ), 
                      go.Surface(z=v[:,:,25],x=np.exp(Kd_mat[:,:,0]), y=Kg_mat[:,:,0], showscale=False), 
                      go.Surface(z=v[:,:,-1],x=np.exp(Kd_mat[:,:,0]), y=Kg_mat[:,:,0], showscale=False), 
                     ])

fig.update_layout(title='V(Kd, Kg), for given Y',
                  width=1000, height=800,
                  margin=dict(l=65, r=50, b=10, t=90),

scene = dict(
                    xaxis_title='Total capital',
                    yaxis_title='Green capital / total capital ratio',
                    zaxis_title='Value function'),
                 )

fig.show()
In [3]:
fig = go.Figure(data=[go.Surface(z=i_d[:,:,0], x=np.exp(Kd_mat[:,:,0]), y=Kg_mat[:,:,0], 
                                 showlegend=True,
                                 name="Y=0",
                                 showscale=False ), 
                      go.Surface(z=i_d[:,:,20],x=np.exp(Kd_mat[:,:,0]), y=Kg_mat[:,:,0], showlegend=True,
                                 name="Y=2",showscale=False), 
                      go.Surface(z=i_d[:,:,-1],x=np.exp(Kd_mat[:,:,0]), y=Kg_mat[:,:,0], name="Y=4",
                                 showlegend=True,
                                 showscale=False), 
                     ])

fig.update_layout(title='i_d(K, r), for given Y',
                  width=1000, height=800,
                  margin=dict(l=65, r=50, b=10, t=90),
                  legend=dict(
    yanchor="top",
    y=0.99,
    xanchor="left",
    x=0.01),
scene = dict(
                    xaxis_title='Total capital',
                    yaxis_title='Green capital / total capital ratio',
                    zaxis_title='Dirty investment'),
                 )

fig.show()
In [4]:
fig = go.Figure(data=[go.Surface(z=i_g[:,:,0], x=np.exp(Kd_mat[:,:,0]), y=Kg_mat[:,:,0], 
                                 showlegend=True,
                                 name="Y=0",
                                 showscale=False ), 
                      go.Surface(z=i_g[:,:,20],x=np.exp(Kd_mat[:,:,0]), y=Kg_mat[:,:,0], showlegend=True,
                                 name="Y=2",showscale=False), 
                      go.Surface(z=i_g[:,:,-1],x=np.exp(Kd_mat[:,:,0]), y=Kg_mat[:,:,0], name="Y=4",
                                 showlegend=True,
                                 showscale=False), 
                     ])

fig.update_layout(title='i_g(K, r), for given Y',
                  width=1000, height=800,
                  margin=dict(l=65, r=50, b=10, t=90),
                  legend=dict(
    yanchor="top",
    y=0.99,
    xanchor="left",
    x=0.01),
scene = dict(
                    xaxis_title='Total capital',
                    yaxis_title='Green capital / total capital ratio',
                    zaxis_title='Green investment'),
                 )

fig.show()
In [5]:
plt.plot(np.exp(Kd )* Kg[2], v[:, 2,25], label="Green capital / Total capital = ".format(Kg[2]))
plt.plot(np.exp(Kd) * Kg[25], v[:, 25,25], label="Green capital / Total capital = ".format(Kg[25]))
plt.plot(np.exp(Kd) * Kg[-1], v[:, -1,25], label="Green capital / Total capital = ".format(Kg[-1]))
plt.title("Value function")
plt.xlabel("Green capital")
Out[5]:
Text(0.5, 0, 'Green capital')
In [6]:
plt.plot(np.exp(Kd), i_d[:, 0,20], label="Clean capital / Total capital = {}".format(Kg[0]) )
plt.plot(np.exp(Kd), i_d[:, 25,20],  label="Clean capital / Total capital = {}".format(Kg[25]))
plt.plot(np.exp(Kd), i_d[:, -1,20],  label="Clean capital / Total capital = {}".format(Kg[-1]))
plt.ylim(0)
plt.xlabel("Total capital")
plt.title("Dirty investment, Y = 2")
plt.legend()
Out[6]:
<matplotlib.legend.Legend at 0x7f9dabe98be0>
In [7]:
plt.plot(np.exp(Kd), i_d[:, 0,40], label="Clean capital / Total capital = {}".format(Kg[0]))
plt.plot(np.exp(Kd), i_d[:, 25,40], label="Clean capital / Total capital = {}".format(Kg[25]))
plt.plot(np.exp(Kd), i_d[:, -1,40], label="Clean capital / Total capital = {}".format(Kg[-1]))
plt.xlabel("Total capital")
plt.title("Dirty investment, Y = 4")
plt.ylim(0)
plt.legend()
Out[7]:
<matplotlib.legend.Legend at 0x7f9dabcfd310>
In [8]:
plt.plot(Y, v[25, 0, :], label="Green capital / Total capital = ".format(Kd[0]))
plt.plot(Y, v[25, 25,:], label="Green capital / Total capital = ".format(Kd[25]))
plt.plot(Y, v[25, -1,:], label="Green capital / Total capital = ".format(Kd[-1]))
# plt.ylim(0)
plt.title("Value function")
plt.xlabel("Temperature anomaly")
Out[8]:
Text(0.5, 0, 'Temperature anomaly')
In [9]:
plt.plot(Y, i_d[0, 25, :])
plt.plot(Y, i_d[25, 25,:])
plt.plot(Y, i_d[-1, 25,:])
plt.title("Dirty investment")
plt.xlabel("Temperature anomaly")
plt.ylim(0)
Out[9]:
(0.0, 0.08542750080854612)
In [10]:
# fig 1
plt.figure()
plt.title("value function - total capital")
plt.plot(np.exp(Kd), v[:, 10,10], label="Green capital / Total capital = {:.2f}, $Y = {:.2f}$".format(Kg[10], Y[10]))
plt.xlabel("Total capital")
plt.ylabel("v")
plt.legend()
# plt.savefig(figuredir + "v-Kd.pdf")
plt.show()
In [11]:
# fig 2
plt.figure()
plt.title("value function - green capital")
plt.plot(np.exp(Kd[30]) * Kg, v[30,:,20], label="Total capital = {:d}, $Y = {:.2f}$".format(int(np.exp(Kd[30])), Y[20]))
plt.xlabel("Green capital")
plt.ylabel("v")
plt.legend()
# plt.savefig(figuredir + "v-Kg.pdf")
plt.show()
In [12]:
# fig 3
plt.figure()
plt.title("value function - temperature")
plt.plot(Y, v[20,20,:], label="Total capital = {:d},\n green capital / total capital = {:.2f}".format(int(np.exp(Kd[20])), Kg[20]))
# plt.plot(Y, v[5,5,:], label="$K_d = {:d}, K_g = {:d}$".format(int(Kd[5]), int(Kg[5])))
plt.xlabel("$Y$")
plt.ylabel("v")
plt.legend()
# plt.savefig(figuredir + "v-Y.pdf")
plt.show()
In [13]:
# fig 1
plt.figure()
plt.title("dirty investment - total capital")
plt.plot(np.exp(Kd), i_d[:,50,20], 
         label="green capital / total capital = {:.2f}, $Y = {:.2f}$".format(Kg[50], Y[20]))
plt.xlabel("Total capital")
plt.ylabel("$i_d$")
plt.legend()
# plt.savefig(figuredir + "id-Kd.pdf")
plt.show()
In [14]:
# fig 2
plt.figure()
plt.title("dirty investment - green capital / total capital")
plt.plot(Kg, i_d[20,:,20], label="$K_d = {:d}, Y = {:.2f}$".format(int(np.exp(Kd[20])), Y[20]))
plt.xlabel("green capital / total capital")
plt.ylabel("$i_d$")
plt.legend()
# plt.savefig(figuredir + "id-Kg.pdf")
plt.show()
In [15]:
# fig 3
plt.figure()
plt.title("dirty investment - temperature")
plt.plot(Y, i_d[20,20,:], label="$K_d = {:d}, K_g = {:d}$".format(int(np.exp(Kd[20])), int(Kg[20])))
plt.xlabel("$Y$")
plt.ylabel("$i_d$")
plt.legend()
# plt.savefig(figuredir + "id-Y.pdf")
plt.show()
# print(i_g)
In [16]:
# fig 1
plt.figure()
plt.title("green investment - total capital")
plt.plot(np.exp(Kd), i_g[:,50,20], label="green capital / total = {:.2f}, $Y = {:.2f}$".format(Kg[50], Y[20]))
plt.xlabel("total capital")
plt.ylabel("$i_g$")
plt.legend()
# plt.savefig(figuredir + "ig-Kd.pdf")
plt.show()
In [17]:
# fig 2
plt.figure()
plt.title("green investment - total capital")
plt.plot(Kg, i_g[20,:,20], label="total capital = {:d}, $Y = {:.2f}$".format(int(np.exp(Kd[20])), Y[20]))
plt.xlabel("green capital / total capital")
plt.ylabel("$i_g$")
plt.legend()
# plt.savefig(figuredir + "ig-Kg.pdf")
plt.show()
In [18]:
# fig 3
plt.figure()
plt.title("green investment - temperature")
plt.plot(Y, i_g[20,50,:], 
         label="total capital = {:d}, \n green capital / total capital = {:.2f}".format(int(Kd[20]), Kg[50]))
plt.xlabel("$Y$")
plt.ylabel("$i_g$")
plt.legend()
# plt.savefig(figuredir + "ig-Y.pdf")
plt.show()